%% observation probability
function [p] = op(Lambda, X, return_func)

[M, N] = size(X);

% zostaw tylko stany q
Lambda = bitand(Lambda, 1);

% suma obserwacji w stanach q = 0, 1
c = [];
c(2) = sum(sum(Lambda));
c(1) = M .* N - c(2);

% zakladamy wartosci pikseli [0, 255]
X_max = 255;

% wektor obserwacji w stanie q = 0
a_ind = find(~Lambda);
a = X(a_ind);
% wektor obserwacji w stanie q = 1
b_ind = find(Lambda);
b = X(b_ind);

% histogram obserwacji w stanie q = 0 
[a, a_bin] = histc(a, 0:X_max);
% histogram obserwacji w stanie q = 1 
[b, b_bin] = histc(b, 0:X_max);

if nargin > 2 && return_func == true
    %% zwraca wszystkie wartosci funkcji
   p = [a' ./ c(1); b' ./ c(2)];
   p(p == 0) = 1;
%     pa = [a'; b'];
%     pb = c;
else
    %% zwraca tylko wartosci wystepujae w obrazie (gotowe do wymnozenia)
    p = zeros(M, N);
    p(a_ind) = a(a_bin) ./ c(1);
    p(b_ind) = b(b_bin) ./ c(2);
    p = p([2 : M-1], [2 : N-1]);
%      pa = [a(a_bin); b(b_bin)];
%      pb = c;
end

